Multiomics reveals persistence of obesity-associated adipose tissue immune cell phenotypes during weight loss and subsequent weight regain
Figure 5
Matthew A. Cottam
Setup
library(Seurat)
library(RColorBrewer)
library(ggplot2)
library(dplyr)
library(ggpubr)
library(SingleCellExperiment)
library(DESeq2)
library(Matrix.utils)
library(magrittr)
library(stringr)
library(pheatmap)
library(tibble)
library(purrr)
library(tidyr)
library(DEGreport)
library(ggrepel)
library(reticulate)
library(viridis)
#Data file
data.integrated <- readRDS("../Inputs/IntegratedData.rds")
#Color panels
maccols <- brewer.pal(n=8, name="Blues")[c(-1,-3,-5,-7)]
monocols <- c("#ff8ade","#e324ad")
dccols <- brewer.pal(n=9, name="Greens")[-1]
tcols <- brewer.pal(n=8, name="Reds")[-1]
nkcols <- c("#876149","#6e3f22")
bcols <- brewer.pal(n=4, name="Purples")[-1]
othcols <- c("#71a157","#00c5cc","#a18e25","#b282b3")
strcols <- brewer.pal(n=4, name="Oranges")[-1]
wccols = c("#878787", "#518db6","#94cc73","#e96b53")
cols <- c(maccols,monocols,dccols,tcols,nkcols,bcols,othcols,strcols)
llcols <- c("#4292C6","#ff8ade","#238B45","#EF3B2C","#876149","#9E9AC8","#71a157","#00c5cc","#a18e25","#b282b3","#FD8D3C")
#Calculate cellspm
cellspm <- table(data.integrated$individual_mice, data.integrated$highlevel2)
propcellspm <- prop.table(cellspm, margin=1)
#Convert to a dataframe
cellspm <- as.data.frame(cellspm)
propcellspm <- as.data.frame(propcellspm)
#Change column names
colnames(cellspm) <- c("Mouse","Cluster","Counts")
colnames(propcellspm) <- c("Mouse","Cluster","Frequency")
#Add a column indicating the groups to each table
cellspm$Group <- gsub("_Hashtag.","", cellspm$Mouse)
propcellspm$Group <- gsub("_Hashtag.","", propcellspm$Mouse)
#Normalize "Counts" to "Counts per hundred cells sequenced"
cphc<-list()
cellspm <- cellspm %>% arrange(.,Mouse)
for (i in levels(cellspm$Mouse)){
x <- sum(subset(cellspm, Mouse==i)$Counts)/100
y <- subset(cellspm, Mouse==i)$Counts/x
cphc <- append(cphc, y)
}
cellspm$CountsPerHundred <- as.numeric(cphc)
#Function to calculate standard error for each given variable (cell cluster)
data_summary <- function(data, varname, groupnames){
require(plyr)
summary_func <- function(x, col){
c(mean = mean(x[[col]], na.rm=TRUE),
se = sd(x[[col]], na.rm=TRUE) / sqrt(sum(!is.na(x[[col]]))))
}
data_sum<-ddply(data, groupnames, .fun=summary_func,
varname)
names(data_sum)[names(data_sum) == 'mean'] <- varname
return(data_sum)
}
#Summarize the data
clustercounts <- data_summary(cellspm, varname="CountsPerHundred", groupnames=c("Group","Cluster"))
clusterprops <- data_summary(propcellspm, varname="Frequency", groupnames=c("Group","Cluster"))Figure 5
Figure 5A
fig5cols <- c(rep("#e6e6e6",4),"#ff8ade","#e324ad","#E5F5E0","#C7E9C0","#A1D99B","#74C476","#41AB5D","#238B45","#006D2C","#00441B",rep("#e6e6e6",19))
DimPlot(data.integrated, split.by="orig.ident",group.by="highlevel2") +
NoLegend() +
labs(title="") &
scale_color_manual(values=fig5cols) &
NoAxes()
Figure 5B
Monocytes <- subset(data.integrated, highlevel2=="Classical Monocytes" | highlevel2=="Non-classical Monocytes")
Monocytes$highlevel2 <- factor(Monocytes$highlevel2, levels=c("Classical Monocytes","Non-classical Monocytes"))Idents(Monocytes) <- Monocytes$highlevel2
DefaultAssay(Monocytes) <- "RNA"
VlnPlot(Monocytes,
features=c("Lyz2","Lst1","Plac8","Fcgr1","Ly6c2","Ccr2","Cd14","Ear2","Cx3cr1","Ace"),
stack=T,
flip=T,
fill.by="ident",
cols=c("#ff8ade","#e324ad","#00441B","#E5F5E0","#74C476")) +
labs(x="") +
NoLegend() +
theme(aspect.ratio = 0.25) +
scale_x_discrete(labels=c("Classical","Non-classical"))
Figure 5C
monocellspm <- subset(cellspm,
Cluster=="Classical Monocytes" |
Cluster=="Non-classical Monocytes" )
monocellspm$Group <- factor(monocellspm$Group, levels = c("Lean","Obese","WL","WC"))
monoclustercounts <- subset(clustercounts,
Cluster=="Classical Monocytes" |
Cluster=="Non-classical Monocytes")
monoclustercounts$Group <- factor(monoclustercounts$Group, levels = c("Lean","Obese","WL","WC"))#T tests with bonferroni correction are used here comparing all groups to the reference lean control group.
compare_means(CountsPerHundred ~ Group, data = subset(monocellspm, Cluster=="Classical Monocytes"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(monocellspm, Cluster=="Non-classical Monocytes"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")ggboxplot(d=subset(monocellspm, Cluster=="Classical Monocytes"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank(), strip.text = element_text(face="bold", size=10)) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,20)
ggboxplot(d=subset(monocellspm, Cluster=="Non-classical Monocytes"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank(), strip.text = element_text(face="bold", size=10)) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,5)
#Just for the legend
ggboxplot(d=subset(monocellspm, Cluster=="Non-classical Monocytes"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.x=element_text(face="bold", size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
ylim(0,5)
Figure 5D
#Convert Seurat Object into a SingleCellExperiment object
counts <- data.integrated@assays$RNA@counts #Raw counts
metadata <- data.integrated@meta.data[,c("orig.ident","HTO_maxID","HTO_classification.global","lowlevel2","highlevel2","individual_mice")] #metadata
metadata$individual_mice <- gsub("_","",metadata$individual_mice)
sce <- SingleCellExperiment(assays=list(counts=counts), colData = metadata)
#Identify Groups to aggregate counts from
groups <- colData(sce)[, c("lowlevel2", "individual_mice")]
#Aggregrate across cluster-sample groups
pb <- aggregate.Matrix(t(counts(sce)),
groupings = groups, fun = "sum")
#Create a vector to split samples.
splitf <- sapply(stringr::str_split(rownames(pb),
pattern = "_",
n = 2),
`[`, 1)
#Split matrix into a list containing components for each cluster and transpose
pb <- split.data.frame(pb,
factor(splitf)) %>%
lapply(function(u)
set_colnames(t(u),
stringr::str_extract(rownames(u), "(?<=_)[:alnum:]+")))
#Create vectors containing sample information
clusters <- purrr::set_names(levels(sce$lowlevel2))
n_clusters <- length(clusters)
samples <- purrr::set_names(levels(factor(sce$individual_mice)))
n_samples <- length(samples)
#Collect sample metadata
n_cells <- as.numeric(table(sce$individual_mice))
m <- match(samples, sce$individual_mice)
ei <- data.frame(colData(sce)[m,],
n_cells, row.names=NULL) %>%
select(-"lowlevel2")
#Create a table to see how many cells belong to each grouping
data.frame(table(sce$individual_mice, sce$lowlevel2)) %>% 'colnames<-'(c("Sample","Subcluster","# Cells")) %>% tibble()#Generate sample-level metadata
metadata <- unique(data.frame(metadata$lowlevel2, metadata$individual_mice, metadata$orig.ident)) %>% 'colnames<-'(c("cluster_id","sample_id","group_id"))
#Vector of cluster IDs
clusters <- levels(metadata$cluster_id)
#Subset for the group of interest
cluster_metadata <- metadata[which(metadata$cluster_id == "Monocytes"), ]
rownames(cluster_metadata) <- cluster_metadata$sample_id
counts <- pb$`Monocytes`
cluster_counts <- data.frame(counts[, which(colnames(counts) %in% rownames(cluster_metadata))])
cluster_metadata <- cluster_metadata[match(colnames(cluster_counts),cluster_metadata$sample_id),]
all(rownames(cluster_metadata) == colnames(cluster_counts)) ## [1] TRUE#Create the DESeq2 data object
dds <- DESeqDataSetFromMatrix(countData=round(cluster_counts),
colData = cluster_metadata,
design = ~ group_id)## converting counts to integer mode#Set obese as the reference control. We are using LRT insteaed of Wald, so this doesn't really matter.
dds$group_id <- relevel(dds$group_id, ref="Obese")
rld <- rlog(dds, blind=TRUE)
#Plot PCA
DESeq2::plotPCA(rld, intgroup = "group_id")
#Plot heatmap
rld_mat <- assay(rld)
rld_cor <- cor(rld_mat)
pheatmap(rld_cor, annotation = cluster_metadata[, c("group_id"), drop=F])
#Run LRT testing via DESeq2
dds <- DESeq(dds, test = "LRT", reduced=~1)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testing#plot dispersion
plotDispEsts(dds)
#Define contrasts. LRT compares across all groups, so this isn't necessary either.
contrast <- c("group_id", levels(cluster_metadata$group_id)[1], levels(cluster_metadata$group_id)[2])
#Generate DE results
res <- results(dds,
contrast = contrast,
alpha = 0.05)
res_tbl <- res %>%
data.frame() %>%
rownames_to_column(var="gene") %>%
as_tibble()
#Filter the results, keeping only significant genes
padj_cutoff <- 0.05
sig_res <- dplyr::filter(res_tbl, padj < padj_cutoff) %>%
dplyr::arrange(padj)
#Normalize the counts for all genes
normalized_counts <- counts(dds,
normalized = TRUE)
#Plot the top 20 genes
top20_sig_genes <- sig_res %>%
dplyr::arrange(padj) %>%
dplyr::pull(gene) %>%
head(n=20)
top20_sig_norm <- data.frame(normalized_counts) %>%
rownames_to_column(var = "gene") %>%
dplyr::filter(gene %in% top20_sig_genes)
gathered_top20_sig <- top20_sig_norm %>%
gather(colnames(top20_sig_norm)[2:length(colnames(top20_sig_norm))], key = "samplename", value = "normalized_counts")
gathered_top20_sig <- inner_join(ei[, c("individual_mice", "orig.ident" )], gathered_top20_sig, by = c("individual_mice" = "samplename"))
ggplot(gathered_top20_sig) +
geom_point(aes(x = gene,
y = normalized_counts,
color = orig.ident),
position=position_jitter(w=0.1,h=0)) +
scale_y_log10() +
xlab("Genes") +
ylab("log10 Normalized Counts") +
ggtitle("Top 20 Significant DE Genes") +
theme_bw() +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(plot.title = element_text(hjust = 0.5)) +
scale_color_manual(values=wccols)
#Use the degreport package to group genes by changes in patterns across groups. Normally this is for timecourse data, but works really effectively in our paradigm to look at genes that recover or fail to recover during WL.
cluster_rlog <- rld_mat[sig_res$gene, ]
metadata <- subset(metadata, cluster_id=="Monocytes")
rownames(metadata) <- metadata$sample_id
metadata <- metadata[match(colnames(cluster_rlog),cluster_metadata$sample_id),]
clusters <- degPatterns(cluster_rlog, metadata = metadata, time = "group_id", col=NULL, plot=FALSE)## Working with 2430 genes.## Working with 2418 genes after filtering: minc > 15## Joining, by = "merge"
## Joining, by = "merge"## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.
## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.ggplot(clusters$normalized, aes(x=group_id,y=value, color=group_id)) +
stat_smooth(aes(x=group_id, y=value))+
geom_point(alpha=0.5, size=1, position=position_jitterdodge(dodge.width = 0.9)) +
geom_line(aes(group=genes), alpha=0.1) +
geom_boxplot(alpha=0, outlier.size=0, outlier.shape=NA)+
facet_wrap(~cluster, ncol=4) +
theme_pubr(base_size = 12) +
labs(y="Z-score of gene abundance") +
scale_color_manual(values=wccols)+
theme(axis.title.x=element_blank()) +
NoLegend()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#Gene clusters:
# 1 = Lean-associated genes that recover with WL
# 2 = Obese-associated genes that do not recover with WL
# 3 = Obese-associated genes that recover with WL
# 4 = Obese-associated genes that do not recover with WL
# 5 = Lean-associated genes that do not recover with WL
# 6 = Lean-associated genes that do not recover with WL
# 7 = Obese-associated genes that do not recover with WL
#Regroups clusters by the above simplified patterns into.
Monocyte_DEG_patterns <- clusters$normalized %>% mutate(patterns = ifelse(cluster==1, 1, #Highest in lean, recover with WL
ifelse(cluster==3,2, #Lowest in lean, recover with WL
ifelse(cluster==2 | cluster==4 | cluster==7,4, #lowest in lean, don't recover with WL
ifelse(cluster==5|cluster==6,3,""))))) #highest in lean, don't recover with WL
ggplot(Monocyte_DEG_patterns, aes(x=group_id,y=value, color=group_id)) +
stat_smooth(aes(x=group_id, y=value))+
geom_point(alpha=0.5, size=0.5, position=position_jitterdodge(dodge.width = 0.9)) +
geom_line(aes(group=genes), alpha=0.1) +
geom_boxplot(alpha=0, outlier.size=0, outlier.shape=NA)+
facet_wrap(~patterns, ncol=4) +
theme_pubr(base_size = 12) +
labs(y="Z-score of gene abundance") +
scale_color_manual(values=wccols)+
theme(strip.background = element_blank(), strip.text.x=element_blank(), axis.title.x=element_blank()) +
NoLegend()## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
#To save genes for putting into Metascape
#write.csv(data.frame(Monocyte_DEG_patterns$genes, Monocyte_DEG_patterns$patterns) %>% unique() %>% 'colnames<-'(c("genes","patterns")), file="Monocyte_metascape.csv")#Ran these through metascape for gene ontology
Pattern1 <- read.csv("../Inputs/Fig5/Monocyte_Pattern_A/Enrichment_GO/_FINAL_GO.csv")
Pattern2 <- read.csv("../Inputs/Fig5/Monocyte_Pattern_B/Enrichment_GO/_FINAL_GO.csv")
Pattern3 <- read.csv("../Inputs/Fig5/Monocyte_Pattern_C/Enrichment_GO/_FINAL_GO.csv")
Pattern4 <- read.csv("../Inputs/Fig5/Monocyte_Pattern_D/Enrichment_GO/_FINAL_GO.csv")
#We will plot the top5 pathways for each pattern
top_pathways_1 <- subset(Pattern1,FirstInGroupByEnrichment==1)
top_pathways_1$GO <- factor(top_pathways_1$GO, levels = top_pathways_1$GO[order(top_pathways_1$Log.q.value.)])
top_pathways_1 <- top_pathways_1 %>% distinct(Description,.keep_all = TRUE)
ggplot(top_pathways_1[1:5,],aes(x=-Log.q.value., y=rev(GO), color=Log.q.value., label=stringr::str_wrap(Description,30)))+
geom_point(size=4) +
labs(x="-log10(adj P value)",y="") +
theme_pubr(base_size=10) +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
xlim(0,30) +
scale_color_gradient(high="#d16f6f",low="#960000") +
NoLegend() +
geom_text_repel(nudge_x = 1,nudge_y = 0.55, size=3, min.segment.length = 0.001) +
geom_vline(xintercept=-log10(0.01), color="grey",linetype="dashed")
top_pathways_2 <- subset(Pattern2,FirstInGroupByEnrichment==1)
top_pathways_2$GO <- factor(top_pathways_2$GO, levels = top_pathways_2$GO[order(top_pathways_2$Log.q.value.)])
top_pathways_2 <- top_pathways_2 %>% distinct(Description,.keep_all = TRUE)
ggplot(top_pathways_2[1:5,],aes(x=-Log.q.value., y=rev(GO), color=Log.q.value., label=stringr::str_wrap(Description,20)))+
geom_point(size=4) +
labs(x="-log10(adj P value)",y="") +
theme_pubr(base_size=10) +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
xlim(0,10) +
scale_color_gradient(high="#d16f6f",low="#960000") +
NoLegend() +
geom_text_repel(nudge_x = 3,nudge_y = 0.33, size=3, min.segment.length = 0.001) +
geom_vline(xintercept=-log10(0.01), color="grey",linetype="dashed")
top_pathways_3 <- subset(Pattern3,FirstInGroupByEnrichment==1)
top_pathways_3$GO <- factor(top_pathways_3$GO, levels = top_pathways_3$GO[order(top_pathways_3$Log.q.value.)])
top_pathways_3 <- top_pathways_3 %>% distinct(Description,.keep_all = TRUE)
ggplot(top_pathways_3[1:5,],aes(x=-Log.q.value., y=rev(GO), color=Log.q.value., label=stringr::str_wrap(Description,25)))+
geom_point(size=4) +
labs(x="-log10(adj P value)",y="") +
theme_pubr(base_size=10) +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
xlim(0,10) +
scale_color_gradient(high="#d16f6f",low="#960000") +
NoLegend() +
geom_text_repel(nudge_x = -4,nudge_y = 0.6, size=3, min.segment.length = 0.001) +
geom_vline(xintercept=-log10(0.01), color="grey",linetype="dashed")
top_pathways_4 <- subset(Pattern4,FirstInGroupByEnrichment==1)
top_pathways_4$GO <- factor(top_pathways_4$GO, levels = top_pathways_4$GO[order(top_pathways_4$Log.q.value.)])
top_pathways_4 <- top_pathways_4 %>% distinct(Description,.keep_all = TRUE)
ggplot(top_pathways_4[1:5,],aes(x=-Log.q.value., y=rev(GO), color=Log.q.value., label=stringr::str_wrap(Description,20)))+
geom_point(size=4) +
labs(x="-log10(adj P value)",y="") +
theme_pubr(base_size=10) +
theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
xlim(0,30) +
scale_color_gradient(high="#d16f6f",low="#960000") +
NoLegend() +
geom_text_repel(nudge_x = -6,nudge_y = 0.5, size=3, min.segment.length = 0.001) +
geom_vline(xintercept=-log10(0.01), color="grey",linetype="dashed")
Figure 5E
Fig5DCs <- subset(data.integrated, highlevel2=="cDC1" | highlevel2=="Activated cDC1" | highlevel2=="Cycling cDC1" | highlevel2=="cDC2" | highlevel2=="Activated cDC2" | highlevel2=="Cycling cDC2" | highlevel2=="pDCs" | highlevel2=="moDCs")Idents(Fig5DCs) <- Fig5DCs$highlevel2
DefaultAssay(Fig5DCs) <- "RNA"
VlnPlot(Fig5DCs,
features=c("Cst3","Clec9a","Xcr1","Sirpa","Cd209a","Stmn1","Pclaf","Siglech","Ear2"),
stack=T,
flip=T,
fill.by="ident",
cols=c("#E5F5E0","#C7E9C0","#A1D99B","#74C476","#41AB5D","#238B45","#006D2C","#00441B")) +
labs(x="") +
NoLegend() +
theme(aspect.ratio = 0.25) 
Figure 5F
dccomps <- list(c("Lean","WL"),c("Lean","Obese"),c("Obese","WC"))
DCscellspm <- subset(cellspm,
Cluster=="cDC1" |
Cluster=="Activated cDC1" |
Cluster=="Cycling cDC1" |
Cluster=="cDC2" |
Cluster=="Activated cDC2" |
Cluster=="Cycling cDC2" |
Cluster=="moDCs" |
Cluster=="pDCs")
DCscellspm$Group <- factor(DCscellspm$Group, levels = c("Lean","Obese","WL","WC"))
DCsclustercounts <- subset(clustercounts,
Cluster=="cDC1" |
Cluster=="Activated cDC1" |
Cluster=="Cycling cDC1" |
Cluster=="cDC2" |
Cluster=="Activated cDC2" |
Cluster=="Cycling cDC2" |
Cluster=="moDCs" |
Cluster=="pDCs")
DCsclustercounts$Group <- factor(DCsclustercounts$Group, levels = c("Lean","Obese","WL","WC"))#T-test with bonferroni correction for multiple comparisons was used to compare all groups to the Lean reference control.
compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="cDC1"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="Activated cDC1"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="Cycling cDC1"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="cDC2"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="Activated cDC2"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="Cycling cDC2"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="moDCs"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")compare_means(CountsPerHundred ~ Group, data = subset(DCscellspm, Cluster=="pDCs"), method="t.test", p.adjust.method = "bonferroni", ref.group="Lean")ggboxplot(d=subset(DCscellspm, Cluster=="cDC1"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,10)
ggboxplot(d=subset(DCscellspm, Cluster=="Activated cDC1"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,2)
ggboxplot(d=subset(DCscellspm, Cluster=="Cycling cDC1"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,2)
ggboxplot(d=subset(DCscellspm, Cluster=="cDC2"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,6)
ggboxplot(d=subset(DCscellspm, Cluster=="Activated cDC2"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,3)
ggboxplot(d=subset(DCscellspm, Cluster=="Cycling cDC2"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,1)
ggboxplot(d=subset(DCscellspm, Cluster=="moDCs"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,8)
ggboxplot(d=subset(DCscellspm, Cluster=="pDCs"),
x="Group",
y="CountsPerHundred",
fill="Group",
facet.by="Cluster",
add = "mean_se") +
theme_classic() +
scale_fill_manual(values=wccols) +
theme(axis.text.y=element_text(face="bold", size=10), strip.text = element_text(face="bold",size=10)) +
labs(x="", y="") +
theme(axis.text.x = element_blank(), axis.ticks = element_blank()) +
facet_wrap(~Cluster) +
NoLegend() +
ylim(0,1)
Figure 5G
Fig5DCs <- subset(data.integrated, highlevel2=="cDC1" | highlevel2=="Activated cDC1" | highlevel2=="Cycling cDC1" | highlevel2=="cDC2" | highlevel2=="Activated cDC2" | highlevel2=="Cycling cDC2" | highlevel2=="pDCs" | highlevel2=="moDCs")Idents(Fig5DCs) <- Fig5DCs$highlevel2
DefaultAssay(Fig5DCs) <- "RNA"
VlnPlot(Fig5DCs,
features=c("Ccr7","Mreg","Fscn1","Il15ra","Il12b","Fabp5","Ldha"),
stack=T,
flip=T,
fill.by="ident",
cols=c("#E5F5E0","#C7E9C0","#A1D99B","#74C476","#41AB5D","#238B45","#006D2C","#00441B")) +
labs(x="") +
NoLegend() +
theme(aspect.ratio = 0.3) 
Figure 5H
For plotting rna velocity, I have a conda env containing scvelo (v0.2.3) and scanpy (v1.7.2). Knitting this whole file takes >~8gb of RAM. I’ve set eval=FALSE in the chunk options so that I can knit this on my local machine quickly. For running this fresh, just change to eval=TRUE.
#Prep environment for reticulate (i.e. we are running Python through R)
use_condaenv("r-velocity", required = TRUE)
conda_list()
scv <- import("scvelo")
scanpy <- import("scanpy")
matplotlib <- import("matplotlib")
plt <- import("matplotlib.pyplot", as = "plt")
#Subset for cDC clusters
Idents(data.integrated) <- data.integrated$lowlevel2
DCvelo <- subset(data.integrated, idents=c("Dendritic Cells"))
Idents(DCvelo) <- DCvelo$highlevel2
DCvelo <- subset(DCvelo, idents=c("cDC1","Activated cDC1","Cycling cDC1","cDC2","Activated cDC2","Cycling cDC2"))
#Create an adata object. Making this from pieces was the only way I could get it to work (Seurat -> adata conversion is difficult).
spliced = DCvelo@assays$spliced@counts
unspliced = DCvelo@assays$unspliced@counts
row.num <- which(rownames(spliced) %in% rownames(DCvelo@assays$RNA@counts))
spliced <- spliced[c(row.num),]
unspliced <- unspliced[c(row.num),]
ad <- import("anndata", convert=FALSE)
orig.ident <- DCvelo$orig.ident
HTO_maxID <- DCvelo$HTO_maxID
lowlevel2 <- DCvelo$lowlevel2
highlevel2 <- DCvelo$highlevel2
clusters <- DCvelo$highlevel2
dfobs <- data.frame(orig.ident, HTO_maxID, lowlevel2, highlevel2, clusters)
rownames(dfobs) <- names(DCvelo$orig.barcodes)
genes_attributes <- rownames(DCvelo@assays$RNA@counts)
dfvar <- data.frame(genes_attributes)
rownames(dfvar) <- rownames(DCvelo@assays$RNA@counts)
emb <- Embeddings(DCvelo, "umap")
adata_DC <- ad$AnnData(
X=t(DCvelo@assays$RNA@counts),
obs=dfobs,
var=dfvar,
layers=list('spliced'=t(spliced), 'unspliced'=t(unspliced)),
obsm=list('X_umap'=emb))
adata_DC
#Run through the scvelo pipeline and generate a dynamic velocity estimate.
scv$pp$filter_genes(adata_DC)
scv$pp$moments(adata_DC)
scv$tl$recover_dynamics(adata_DC)
scv$tl$velocity(adata_DC, mode='dynamical')
scv$tl$velocity_graph(adata_DC)
#Colors for plot
DC.colors <- c("#E5F5E0","#C7E9C0","#A1D99B","#74C476","#41AB5D","#238B45")#These will open a separate window because they run through reticulate.
#save as "scvelo1.png" - This is already available in the "figures" folder.
scv$pl$velocity_embedding_stream(adata_DC,
basis='umap',
size=100,
legend_fontsize=14,
palette=DC.colors,
smooth=TRUE,
min_mass=0,
alpha=0.5)knitr::include_graphics("figures/scvelo1.png")
Figure 5I
DCs <- subset(data.integrated,lowlevel2=="Dendritic Cells")
FeaturePlot(DCs, features=c("Cd274","Pdcd1lg2","Cd200"), split.by="orig.ident", by.col=T, pt.size=1, order=T) & theme_void() & NoLegend() & scale_color_viridis(option="D", limits=c(0,2)) & ylim(-10,3) & xlim(-7,7)## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 79 rows containing missing values (geom_point).## Warning: Removed 30 rows containing missing values (geom_point).## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 46 rows containing missing values (geom_point).## Warning: Removed 79 rows containing missing values (geom_point).## Warning: Removed 30 rows containing missing values (geom_point).## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 46 rows containing missing values (geom_point).## Warning: Removed 79 rows containing missing values (geom_point).## Warning: Removed 30 rows containing missing values (geom_point).## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 46 rows containing missing values (geom_point).
### for scale bar
FeaturePlot(DCs, features=c("Cd274","Pdcd1lg2","Cd200"), split.by="orig.ident", by.col=T, pt.size=0.5) & scale_color_viridis(option="D", limits=c(0,2), breaks=c(0,1,2)) & ylim(-10,3) & xlim(-7,7) & theme(legend.position="bottom")## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.
## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.## Warning: Removed 79 rows containing missing values (geom_point).## Warning: Removed 30 rows containing missing values (geom_point).## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 46 rows containing missing values (geom_point).## Warning: Removed 79 rows containing missing values (geom_point).## Warning: Removed 30 rows containing missing values (geom_point).## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 46 rows containing missing values (geom_point).## Warning: Removed 79 rows containing missing values (geom_point).## Warning: Removed 30 rows containing missing values (geom_point).## Warning: Removed 64 rows containing missing values (geom_point).## Warning: Removed 46 rows containing missing values (geom_point).
### Just using this for nice labels
VlnPlot(DCs, features=c("Cd274","Pdcd1lg2","Cd200"), stack=T, flip=T, split.by="orig.ident")## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
Session Info
sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] plyr_1.8.6 viridis_0.6.1
## [3] viridisLite_0.4.0 reticulate_1.20
## [5] ggrepel_0.9.1 DEGreport_1.28.0
## [7] tidyr_1.1.3 purrr_0.3.4
## [9] tibble_3.1.2 pheatmap_1.0.12
## [11] stringr_1.4.0 magrittr_2.0.1
## [13] Matrix.utils_0.9.8 Matrix_1.3-4
## [15] DESeq2_1.32.0 SingleCellExperiment_1.14.1
## [17] SummarizedExperiment_1.22.0 Biobase_2.52.0
## [19] GenomicRanges_1.44.0 GenomeInfoDb_1.28.0
## [21] IRanges_2.26.0 S4Vectors_0.30.0
## [23] BiocGenerics_0.38.0 MatrixGenerics_1.4.0
## [25] matrixStats_0.59.0 ggpubr_0.4.0
## [27] dplyr_1.0.7 ggplot2_3.3.5
## [29] RColorBrewer_1.1-2 SeuratObject_4.0.2
## [31] Seurat_4.0.3
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.1 tidyselect_1.1.1
## [3] RSQLite_2.2.7 AnnotationDbi_1.54.1
## [5] htmlwidgets_1.5.3 grid_4.1.0
## [7] BiocParallel_1.26.0 Rtsne_0.15
## [9] munsell_0.5.0 codetools_0.2-18
## [11] ica_1.0-2 future_1.21.0
## [13] miniUI_0.1.1.1 withr_2.4.2
## [15] colorspace_2.0-2 highr_0.9
## [17] knitr_1.33 ROCR_1.0-11
## [19] ggsignif_0.6.2 tensor_1.5
## [21] listenv_0.8.0 labeling_0.4.2
## [23] lasso2_1.2-21.1 GenomeInfoDbData_1.2.6
## [25] mnormt_2.0.2 polyclip_1.10-0
## [27] farver_2.1.0 bit64_4.0.5
## [29] parallelly_1.27.0 vctrs_0.3.8
## [31] generics_0.1.0 xfun_0.24
## [33] doParallel_1.0.16 R6_2.5.0
## [35] clue_0.3-59 locfit_1.5-9.4
## [37] reshape_0.8.8 bitops_1.0-7
## [39] spatstat.utils_2.2-0 cachem_1.0.5
## [41] DelayedArray_0.18.0 assertthat_0.2.1
## [43] promises_1.2.0.1 scales_1.1.1
## [45] gtable_0.3.0 Cairo_1.5-12.2
## [47] globals_0.14.0 goftest_1.2-2
## [49] rlang_0.4.11 genefilter_1.74.0
## [51] GlobalOptions_0.1.2 splines_4.1.0
## [53] rstatix_0.7.0 lazyeval_0.2.2
## [55] spatstat.geom_2.2-2 broom_0.7.7
## [57] yaml_2.2.1 reshape2_1.4.4
## [59] abind_1.4-5 backports_1.2.1
## [61] httpuv_1.6.1 tools_4.1.0
## [63] psych_2.1.6 logging_0.10-108
## [65] ellipsis_0.3.2 spatstat.core_2.3-0
## [67] jquerylib_0.1.4 ggdendro_0.1.22
## [69] ggridges_0.5.3 Rcpp_1.0.7
## [71] zlibbioc_1.38.0 RCurl_1.98-1.3
## [73] rpart_4.1-15 deldir_0.2-10
## [75] GetoptLong_1.0.5 pbapply_1.4-3
## [77] cowplot_1.1.1 zoo_1.8-9
## [79] grr_0.9.5 haven_2.4.1
## [81] cluster_2.1.2 data.table_1.14.0
## [83] scattermore_0.7 openxlsx_4.2.4
## [85] circlize_0.4.13 lmtest_0.9-38
## [87] RANN_2.6.1 tmvnsim_1.0-2
## [89] fitdistrplus_1.1-5 hms_1.1.0
## [91] patchwork_1.1.1 mime_0.11
## [93] evaluate_0.14 xtable_1.8-4
## [95] XML_3.99-0.6 rio_0.5.27
## [97] readxl_1.3.1 shape_1.4.6
## [99] gridExtra_2.3 compiler_4.1.0
## [101] KernSmooth_2.23-20 crayon_1.4.1
## [103] htmltools_0.5.1.1 mgcv_1.8-36
## [105] later_1.2.0 geneplotter_1.70.0
## [107] DBI_1.1.1 ComplexHeatmap_2.8.0
## [109] MASS_7.3-54 car_3.0-10
## [111] igraph_1.2.6 forcats_0.5.1
## [113] pkgconfig_2.0.3 foreign_0.8-81
## [115] plotly_4.9.4.1 spatstat.sparse_2.0-0
## [117] foreach_1.5.1 annotate_1.70.0
## [119] bslib_0.2.5.1 XVector_0.32.0
## [121] digest_0.6.27 ConsensusClusterPlus_1.56.0
## [123] sctransform_0.3.2 RcppAnnoy_0.0.18
## [125] spatstat.data_2.1-0 Biostrings_2.60.1
## [127] rmarkdown_2.9 cellranger_1.1.0
## [129] leiden_0.3.8 edgeR_3.34.0
## [131] uwot_0.1.10 curl_4.3.2
## [133] shiny_1.6.0 rjson_0.2.20
## [135] lifecycle_1.0.0 nlme_3.1-152
## [137] jsonlite_1.7.2 carData_3.0-4
## [139] limma_3.48.0 fansi_0.5.0
## [141] pillar_1.6.1 lattice_0.20-44
## [143] Nozzle.R1_1.1-1 KEGGREST_1.32.0
## [145] fastmap_1.1.0 httr_1.4.2
## [147] survival_3.2-11 glue_1.4.2
## [149] zip_2.2.0 iterators_1.0.13
## [151] png_0.1-7 bit_4.0.4
## [153] stringi_1.7.3 sass_0.4.0
## [155] blob_1.2.1 memoise_2.0.0
## [157] irlba_2.3.3 future.apply_1.7.0